智汇华云 | AIOps之动态阈值—SARIMA模型详解（二）
近年来，IT运维人工智能（AIOps）已成为了应对IT系统与日俱增的复杂性的很好的解决方案。AIOps基于大数据、数据分析和机器学习来提供洞察力，并为管理现代基础设施和软件所需的任务提供更高水平的自动化（不依赖于人类操作员）。
因此，AIOps具有巨大的价值。展望未来，AIOps将在IT团队提高效率方面发挥关键作用。它还会使应用复杂的下一代技术成为可能，而且那些技术的复杂性是传统解决方案无法胜任的。
华云数据“智汇华云”专栏将为您奉上“AIOps之动态阈值—SARIMA模型详解”系列文章，本期将为您讲解其第二部分。
本期讲解人
张歆纬
运维开发组研发工程师
在自回归模型中，我们使用过去变量的线性组合来预测。
自回归表示这是对于自身变量的回归。
p阶AR模型可以写成：
y_t=c+ϕ_1 y_(t-1)+ϕ_2 y_(t-2)+⋯+ϕ_p y_(t-p)+ε_t
这里ε_t是白噪声，我们把这个叫做AR(p)模型，p阶自回归模型。
下图展示了AR(1)模型和AR(2)模型：
对于AR(1)模型：
当 ϕ_1=0，yt就是白噪声
当 ϕ_1=1并且c=0，yt就是随机漫步
当 ϕ_1=1并且c≠0，yt就是带偏移量的随机漫步
当ϕ_1<0，yt趋向于在正负值之间振动
我们通常会限制AR模型只用于平稳的数据，所以我们对参数有一些限制：
对于AR(1)模型：-1<ϕ_1<1
对于AR(2)模型：-1<ϕ_2<1，ϕ_1+ϕ_2<1，ϕ_2-ϕ_1<1
对于p>2，参数限制就非常复杂，我们可以用python的包来搞定。
不像AR模型中使用过去的预测变量，MA模型使用过去的预测误差。
y_t=c+ε_t+θ_1 ε_(t-1)+θ_2 ε_(t-2)+⋯+θ_q ε_(t-q)
ε_t是白噪声。我们把这个叫做MA(q)模型，q阶移动平均模型。
下图展示了MA(1)模型和MA(2)模型：
我们可以把任意平稳的AR(p)模型写成MA(∞)模型。比如，我们可以把AR(1)模型写成：
y_t=ϕ_1 y_(t-1)+ε_t
=ϕ_1 (ϕ_1 y_(t-2)+ε_(t-1))+ε_t
=ϕ_1^2 y_(t-2)+ϕ_1 ε_(t-1)+ε_t
=ϕ_1^3 y_(t-3)+ϕ_1^2 ε_(t-2)+ϕ_1 ε_(t-1)+ε_t
有-1<ϕ_1<1，当k变大时，ϕ_1^k会变小。所以我们最终会得到：
y_t=ε_t+ϕ_1 ε_(t-1)+ϕ_1^2 ε_(t-2)+ϕ_1^3 ε_(t-3)+⋯
这是一个MA(∞)模型。
如果我们给MA模型加一些限制，我们可以称MA模型是可逆的，我们可以把任意MA(q)模型写成AR(∞)模型。
可逆性限制和平稳性限制类似：
对于MA(1)模型：-1<θ_1<1
对于MA(2)模型：-1<θ_2<1，θ_2+θ_1>-1，θ_1-θ_2<1
对于q>2，参数限制就非常复杂，我们可以用python的包来搞定。
ARIMA模型
Auto Regressive Integrated Moving Average
如果我们组合AR和MA模型并差分，我们可以得到ARIMA模型。
模型可以写成：
y_t^'=c+ϕ_1 y_(t-1)^'+⋯+ϕ_p y_(t-p)^'+θ_1 ε_(t-1)+⋯+θ_q ε_(t-q)+ε_t
y_t^'是差分过的序列，右侧的预测器包含延迟yt和延迟误差。我们叫这个ARIMA(p,d,q)模型：
p：自回归阶数
d：差分次数
q：移动平均阶数
有一些特殊的ARIMA模型如下：
白噪声：ARIMA(0,0,0)
随机漫步：ARIMA(0,1,0)
带偏移量的随机漫步：ARIMA(0,1,0)带常数
自回归：ARIMA(p,0,0)
移动平均：ARIMA(0,0,q)
用后移符号，我们可以把ARIMA模型写成：
常数c在长期预测中十分重要：
如果c=0并且d=0，长期预测值会趋向于0
如果c=0并且d=1，长期预测值会趋向于非零常数
如果c=0并且d=2，长期预测值会变成一条直线
如果c≠0并且d=0，长期预测值会趋向于数据的平均值
如果c≠0并且d=1，长期预测值会变成一条直线
如果c≠0并且d=2，长期预测值会变成二次抛物线
自相关系数测量了yt和yt-k的关系。如果yt和yt-1相关，那么yt-1和yt-2肯定也相关。但这样的话，yt和yt-2可能也相关，仅仅只因为他们都跟yt-1相关，而不是因为yt-2中有新的信息可以用于预测yt。
为了解决这个问题，我们可以使用偏自相关系数。这是在移除延迟1,2,3,…,k-1的影响后，测量yt和yt-k之间的关系。
如果差分过后的ACF和PACF图满足以下形式，数据可能是ARIMA(p,d,0)模型：
ACF是指数衰减或者正弦式的
在PACF中，在延迟p的地方有一个明显的突刺，但后面没有
如果差分过后的ACF和PACF图满足以下形式，数据可能是ARIMA(0,d,q)模型：
PACF是指数衰减或者正弦式的
在ACF中，在延迟q的地方有一个明显的突刺，但后面没有
估算模型的时候，我们使用最大似然估计。已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，参数估计就是通过若干次试验，观察其结果，利用结果推出参数的大概值。
对于ARIMA模型，MLE通过最小化
来获得。
对于给定的p,d,q组合，我们可以用python最大化log likelihood来找到合适的p,d,q。
赤池信息准则(AIC)在选取参数时非常有用，可以写成：
AIC=-2 log(L)+2(p+q+k+1)
其中L是数据的likelihood，如果c=0，k=0；如果c≠0，k=1。
修正赤池信息准则(AICc)可以写成：
AICc=AIC+(2(p+q+k+1)(p+q+k+2))/(T-p-q-k-2)
贝叶斯信息准则(BIC)可以写成：
BIC=AIC+[log(T)-2](p+q+k+1)
最小化AIC，AICc或者BIC可以得到较优模型，我们偏向于选择AIC。
pmdarima是一个python解决ARIMA和SARIMA模型的包，主要使用了Hyndman-Khandakar算法的变形，组合了单位根检验，最小化AICc和MLE。
用于自动化ARIMA模型拟合的Hyndman-Khandakar算法：
1、重复使用KPSS检测决定差分次数
2、差分后最小化AICc来选取p和q的值，这种算法使用了阶梯式搜索来遍历模型空间，而不是考虑所有p和q的组合
a、拟合四个初始模型：
ARIMA(0,d,0)
ARIMA(2,d,2)
ARIMA(1,d,0)
ARIMA(0,d,1)
常数项会被考虑进去除非d=2。如果d≤1，拟合额外的一个模型：
ARIMA(0,d,0)没有常数项
b、在步骤a中最优的模型（最小的AICc值）会被设置为当前模型
c、微调当前模型：
对p或/和q±1
加入/去除常数项c
新的最优模型变成当前模型
d、重复步骤c直到没有更小的AICc
SARIMA模型
Seasonal Auto Regressive Integrated Moving Average
ARIMA模型的缺陷在于没有考虑周期性，加入周期项可以得到SARIMA模型：
ARIMA (p,d,q) (P,D,Q)m
非周期性部分 周期性部分
m是每年的观测数量。P,D,Q作为周期性参数，p,d,q作为非周期性参数。
模型的周期性部分和非周期性部分很相似，但包括了周期后移。
比如，ARIMA(1,1,1)(1,1,1)4对于季度数据(m=4)可以写成：
(1-ϕ_1 B)(1-Φ_1 B^4 )(1-B)(1-B^4 ) y_t=(1+θ_1 B)(1+Θ_1 B^4 ) ε_t
PACF和ACF图中的周期性延迟可以看出AR模型或者MA模型的周期性部分。
比如，SARIMA(0,0,0)(0,0,1)12模型会有以下特性：
ACF中延迟12有突刺，但没有其他的明显突刺
PACF的周期性延迟有指数衰减，比如在延迟12，24，36的地方
相似的，SARIMA(0,0,0)(1,0,0)12模型会有以下特性：
ACF的周期性延迟有指数衰减
PACF中延迟12有突刺
另外，根据简约性原则parsimony principle，p+d+q+P+D+Q≤6为佳。
下面的例子可以很好的解释模型拟合的过程：
例子：欧洲季度零售指数
这个例子是欧洲零售指数从1996到2011年的数据，我们把它套进SARIMA模型进行预测。
这组数据明显是不平稳的，并有一些周期性，所以我们先进行周期性差分，如下图：
这看起来还是不平稳，我们再进行一次差分，如下图：
ACF图中延迟1的明显突刺说明有个非周期性的MA(1)部分，ACF图中延迟4的明显突刺说明有个周期性MA(1)的部分。所以，我们从SARIMA(0,1,1)(0,1,1)4模型开始，得到拟合模型的残差，如下图：
ACF和PACF都在延迟2有明显突刺，延迟3的突刺也不小，所以模型应该还有额外的非周期性部分。SARIMA(0,1,2)(0,1,1)4模型的AICc是74.36，SARIMA(0,1,3)(0,1,1)4模型的AICc是68.53。其他的AR参数都没有更小的AICc值。所以，我们选择SARIMA(0,1,3)(0,1,1)4，画出该模型的残差：
所有突刺都在合理范围内，残差值看起来像白噪声了。Ljung-Box测试也显示残差没有自相关性了。
然后，我们就可以用该模型进行预测了：
图中显示了预测值以及80%和95%的置信区间。
我们已经了解了SARIMA模型，并可以对时间序列数据进行预测了。对于动态阈值，我们首先获取历史数据，对数据进行处理，需要对缺失数据进行一些填充。然后我们进行SARIMA模型拟合，得出最优模型之后，对未来指标走势进行预测，通过95%的置信区间生成阈值区间，如果指标超出这个区间，我们认为指标异常，对用户进行告警。每天我们都会重复以上操作，让模型拟合更加准确，从而使动态阈值功能日趋完善。
相关阅读
智汇华云 | AIOps之动态阈值—SARIMA模型详解（一）
智汇华云 | ping不通 找房东——从一个ping包来解析网络转发原理
点击阅读原文，了解华云数据更多信息